Quantum-to-classical crossover for Andreev billiards in a magnetic field 
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We extend the existing quasiclassical theory for the superconducting proximity effect in a chaotic 
quantum dot, to include a time-reversal-symmetry breaking magnetic field. Random-matrix theory 
(RMT) breaks down once the Ehrenfest time te becomes longer than the mean time td between 
Andreev reflections. As a consequence, the critical field at which the excitation gap closes drops 
below the RMT prediction as te/td is increased. Our quasiclassical results are supported by 
comparison with a fully quantum mechanical simulation of a stroboscopic model (the Andreev 
kicked rotator). 
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I. INTRODUCTION 

When a quantum dot is coupled to a superconductor 
via a point contact, the conversion of electron to hole 
excitations by Andreev reflection governs the low-energy 
spectrum. The density of states of such an Andreev bil- 
liard was calculated using random-matrix theory (RMT) 
. If the classical dynamics in the isolated quantum dot 
is chaotic, a gap opens up in the spectrum. The exci- 
tation gap -Egap is of the order of the Thouless energy 
H/td, with td the average time between Andreev reflec- 
tions. Although chaoticity of the dynamics is essential 
for the gap to open, the size of the gap in RMT is in- 
dependent of the Lyapunov exponent A of the chaotic 
dynamics. 

If the size L of the quantum dot is much larger than 
the Fermi wavelength Xp, a competing timescale te — 
A -1 In (L/Xp) appears, the Ehrenfest time, which causes 
the breakdown of RMT Q . The gap becomes dependent 
on the Lyapunov exponent and for te 3> td vanishes as 
_Eg a p ~ Tl/te- The Ehrenfest time dependence of the gap 
has been investigated in several works 0, 0, IE IE 00 
For a recent review, see Ref. [Toj]. 

A magnetic field breaks time-reversal symmetry, 
thereby reducing -E gap - At a critical field B c the gap 
closes. This was calculated using RMT in Ref. [ll], but 
the effect of a finite Ehrenfest time was not studied be- 
fore. Here we extend the zero-field theory of Silvestrov 
et al. 5] to non-zero magnetic field. It is a quasiclassical 
theory, which relates the excitation spectrum to the clas- 
sical dynamics in the billiard. The entire phase space is 
divided into two parts, depending on the time T between 
Andreev reflections. Times T < te are quantized by 
identifying the adiabatic invariant, while times T > te 
are quantized by an effective RMT with r^-dependent 
parameters. 

There exists an alternative approach to quantization 
of the Andreev billiard, due to Vavilov and Larkin 
which might also be extended to non-zero magnetic field. 
In zero magnetic field the two models have been shown 
to give similar results [H|, so we restrict ourselves here 
to the approach of Ref. 




FIG. 1: Classical trajectory in an Andreev billiard. Particles 
are deflected by the potential V = [(r/L) 2 — ll Vo for r < L, 
V = [-4(r/L) 2 + W(r/L) - 6] V for r > L, with r 2 = x 2 + 
y 2 (the dotted lines are equipotentials) . At the insulating 
boundaries (solid lines) there is specular reflection, while the 
particles are Andreev reflected at the superconductor (y = 0, 
dashed line). Shown is the trajectory of an electron at the 
Fermi level (E = 0), for B = and E F = 0.84 eV . The 
Andreev reflected hole will retrace this path. 



The outline of the paper is as follows. We start by de- 
scribing the adiabatic levels in Sec. [HI followed by the ef- 
fective RMT in Sec. lIIII In Scc. lIVI we compare our quasi- 
classical theory with fully quantum mechanical computer 
simulations. We conclude in SecM 



II. ADIABATIC QUANTIZATION 

We generalize the theory of adiabatic quantization of 
the Andreev billiard of Ref. to include the effect of 
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FIG. 2: Andreev reflection at a NS boundary (dashed line) of 
an electron to a hole. The left panel shows the case of perfect 
retroreflection (zero excitation energy E and zero magnetic 
field B). The middle and right panels show that the hole 
does not precisely retrace the path of the electron if E or B 
are non-zero. 
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FIG. 3: (Color online) Poincare map for the Andreev billiard 
of Fig. Each dot marks the position x and tangential ve- 
locity v x of an electron at the NS boundary. Subsequent dots 
are obtained by following the electron trajectory for E, B — > 
at fixed ratio B/E = |^/m/Voi 2 e 3 . The inset shows the full 
surface of section of the Andreev billiard, while the main plot 
is an enlargement of the central region. The drift is along 
closed contours defined by /C = constant [see Eq. Q]. The 
value of the adiabatic invariant /C (in units of \J mL 2 /eVo) is 
indicated for several contours. All contours are closed loops, 
but for some contours the opening of the loop is not visible 
in the figure. 



a magnetic field. An example of the geometry of such a 
billiard is sketched in Fig.^ The normal metal lies in 
the x-y plane and the boundary with the superconductor 
(NS boundary) is at y = 0. The classical mechanics of 
electrons and holes in such an Andreev billiard has been 
analyzed in Refs. [l2L fl3l IT1 |. We first summarize the 
results we need, then proceed to the identification of the 
adiabatic invariant, and finally present its quantization. 




FIG. 4: Directed area for a classical trajectory, consisting of 
the area enclosed by the trajectory after joining begin and end 
points along the NS boundary (dashed line). Different parts 
of the enclosed area have different signs because the boundary 
is circulated in a different direction. 

A. Classical mechanics 

The classical equation of motion 

f(() = --rxB + -V^(r) (1) 

m m 

is the same for the electron and the hole because both 
charge e and mass m change sign. The vector B is the 
uniform magnetic field in the z-direction and V(r) is the 
electrostatic potential in the plane of the billiard. The 
dots on r = (x, y) denote time derivatives. We follow 
the classical trajectory of an electron starting at the NS 
boundary position (x,0) with velocity (v x ,v y ). The elec- 
tron is at an excitation energy E counted from the Fermi 
level. After a time T the electron returns to the super- 
conductor and is retroreflected as a hole. Retroreflection 
means that v x — > —v x . The y-component v y of the veloc- 
ity also changes sign, but in addition it is slightly reduced 
in magnitude, Vy — > Vy — AE/m, so that an electron at 
an energy E above the Fermi level becomes a hole at an 
energy E below the Fermi level. 

This refraction is one reason why the hole does not 
precisely retrace the path of the electron. A second rea- 
son is that a non-zero B will cause the hole trajectory to 
bend in the direction opposite to the electron trajectory 
(because the velocity has changed sign) , see Fig. [21 It 
follows that if either E or B are non-zero, the hole will 
return to the NS boundary at a slightly different position 
and with a slightly different velocity. The resulting drift 
of the quasi-periodic motion is most easily visualized in 
a Poincare surface of section, see Fig.|3J Each dot marks 
the position x and tangential velocity v x of an electron 
leaving the NS boundary. At non-zero £ or B, subse- 
quent dots are slightly displaced, tracing out a contour 
in the (x, v x ) plane. In the limit E, B — > 0, the shape 
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of these contours is determined by the adiabatic invari- 
ant of the classical dynamics. In Ref. [f| it was shown 
that the contours in the Poincare surface of section are 
isochronous for B = 0. This means that they are given 
by T(x,v x ) — const, with T(x,v x ) the time it takes an 
electron at the Fermi level to return to the NS boundary, 
as a function of the starting point (x,v x ) on the bound- 
ary. In other words, for B — the time between Andreev 
reflections is an adiabatic invariant in the limit E — > 0. 



B. Adiabatic invariant 

We generalize the construction of the adiabatic invari- 
ant of Ref. H| to B 7^ 0. We start from the Poincare 
invariant 



At) 



dr 



(2) 



c(t) 



over a closed contour C(t) in phase space that moves 
according to the classical equations of motion. The con- 
tour extends over two sheets of phase space, joined at 
the NS interface. In the electron sheet the canonical mo- 
mentum is p + = rnv + — eA, while in the hole sheet it 
is p = — mv_ + eA. Both the velocity v±, given in 
absolute value by |v±| = {2/mfl 2 [E F ±E + eF(r)] 1 / 2 
and directed along the motion, as well as the vector po- 
tential A = \Bz x r are functions of the position r on 
the contour, determined, respectively, by the energy E 
and the magnetic field B. (Since the contour is closed, 
the Poincare invariant is properly gauge invariant.) 

Quite generally, dl/dt — 0, meaning that I is a con- 
stant of the motion For E = B = we take C(0) to 
be the self-retracing orbit from electron to hole and back 
to electron. It is obviously time-independent, with 2 = 
(because the contributions from electron and hole sheet 
cancel). For E or B non-zero, we construct C(0) from 
the same closed trajectory in real space, but now with 
p± (r) and A{r) calculated at the given values of E and 
B. Consequently, this contour C(t) will drift in phase 
space, preserving I(t) = 1(0). The Poincare invariant 
is of interest because it is closely related to the action 
integral 



dr. 



(3) 



The action integral is defined as an integral along the pe- 
riodic electron-hole orbit O e h followed by electrons and 
holes at E, B = 0. To every point (x, v x ) in the Poincare 
surface of section corresponds an orbit O e h and hence an 
action integral I(x,v x ). We compare the contour C[t) 
and the trajectory O e h intersecting the Poincare surface 
of section at the same point (x,v x ). At t — they co- 
incide and for sufficiently slow drifts they stay close and 
therefore the action integral I = 1(0) + 0(t 2 ) is an adi- 
abatic invariant of the motion in the Poincare surface of 
section |l5j . 



It remains to determine the adiabatic invariant / in 
terms of E and B and the chosen trajectory C(0). To 
linear order in E, B we find 



I = 2EJC, JC = T - eAB/E, 



(4) 



with A = | §{r x dr) ■ z the directed area (see Fig. 
enclosed by the electron trajectory and the NS bound- 
ary. Both the time T and the area A are to be evaluated 
at E = B = 0. Because E is a constant of the motion, 
adiabatic invariance of I implies that K, = I/2E is an 
adiabatic invariant. At zero field this adiabatic invariant 
is simply the time T between Andreev reflections. At 
non-zero field the invariant time contains also an elec- 
tromagnetic contribution —eAB/E, proportional to the 
enclosed flux. 

Fig. shows that, indeed, the drift in the Poincare 
surface of section is along contours C>c of constant JC. In 
contrast to the zero-field case, the invariant contours in 
the surface of section are now no longer energy indepen- 
dent. This will have consequences for the quantization, 
as we describe next. 



C. Quantization 

The two invariants E and K. define a two-dimensional 
torus in four-dimensional phase space. The two topo- 
logically independent closed contours on this torus are 
formed by the periodic electron-hole orbit O e h and the 
contour Ck. in the Poincare surface of section. The area 
they enclose is quantized following the prescription of 
Einstein-Brillouin-Keller 16, 17], 

f p dr = 2irh(m + l/2), m = 0, 1, 2, . . . (5a) 

f p x dx = 2nh(n + l/2), n= 0,1,2,... (5b) 

JCk. 

The action integral l|5a|l can be evaluated explicitly, lead- 
ing to 



EJC = nh(m+ 1/2). 



(6) 



The second quantization condition l|5bj) gives a second 
relation between E and /C, so that one can eliminate 
/C and obtain a ladder of levels E mn . For B = the 
quantization condition (|5b|) is independent of E, so one 
obtains separately a quantized time T n and quantized en- 
ergy E mn = (m + l/2)Trh/T n . For B ^ both fC mn and 
E m n depend on the sets of integers m, n. 



D. Lowest adiabatic level 

The value E 00 of the lowest adiabatic level follows from 
the pair of quantization conditions (J5J with m = n = 0. 
To determine this value we need to determine the area 
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FIG. 5: Illustration of a bunch of trajectories within a sin- 
gle scattering band in the billiard defined in Fig. All tra- 
jectories in this figure have starting conditions in the band 
containing the contour with K, — 11 of Fig.[3] Both T and A 
vary only slightly from one trajectory to the other, so that the 
whole band can be characterized by a single T and A, being 
the average of T and A over the scattering band. 



O(ZC) = § c ^p x dx enclosed by contours of constant /C, in 
the limit of large K.. 

In Ref. Q the area 0{JC) was determined in the case 
-6 = 0, when K = T and the contours are isochronous. 
It was found that 



0(T) < O exp(-AT), 



(7) 



with A the Lyapunov exponent of the normal billiard 
without superconductor and Oq a characteristic area that 
depends on the angular distribution of the beam of elec- 
trons entering the billiard (width L) from the narrow con- 
tact to the superconductor (width W). For a collimated 
beam having a spread of velocities \v x /vp\ < W/L one 
has Oq = Nh. For a non-collimated beam Oq = NhW/L. 
The integer N is the number of scattering channels con- 
necting the billiard to the superconductor. The quanti- 
zation requirement 0(T) > nh gives the lowest adiabatic 
level in zero magnetic field ||, 



E Q0 {B = 0) 



nh 



T E = - In (OoM) • (8) 



The Ehrenfest time te corresponds to a contour that 
encloses an area nh. 

In order to generalize Eq. J7J to B ^ 0, we discuss 
the concept of scattering bands, introduced in Ref. 
for a normal billiard (where they were called transmis- 
sion and reflection bands) . Scattering bands are ordered 
phase space structures that appear in open systems, even 



if their closed counterparts are fully chaotic. These struc- 
tures are characterized by regions in which the functions 
T(x,v x ) and A(x,v x ) vary slowly almost everywhere. 
Hence, they contain orbits of almost constant return time 
and directed area, that is, orbits returning by bunches. 
One such bunch is depicted in Fig. The scattering 
bands are bounded by contours of diverging T(x, v x ) and 
A(x,v x ). The divergence is very slow (oc 1/lne, with e 
the distance from the contour [J), so the mean return 
time T and mean directed area A in a scattering band 
remain finite and well defined [Tflj . 

The area Oband of a band depends on T as |l8[ 



O b and(T)~O exp(-AT). 



(9) 



Since an isochronous contour must lie within a single 
scattering band, Eq. Q follows from Eq. |J5J and from 
the fact that the distribution of return times is sharply 
peaked around the mean T . Because contours of constant 
JC = T — eAB I E must also lie within a single scattering 
band, the area 0{1C) is bounded by the same function 
Oband (T). We conclude that within a given scattering 
band the largest contour of constant T and the largest 
contour of constant K, each have approximately the same 
area as the band itself, 



0{T),0(K.) < O b and(T) ~ O exp(-AT). 



(10) 



We are now ready to determine the magnetic field de- 
pendence of the lowest adiabatic level Eqq{B). The cor- 
responding contour C/c lies in a band characterized by a 
mean return time T = A -1 In {Oq/ttK), according to Eqs. 
(|5bjl and GDJ. This is the same Ehrenfest time as Eq. 
(JSJl for B — (assuming that the orbital effect of the 
magnetic field does not modify A) . The energy of the 
lowest adiabatic level Eqo is determined by the quantiza- 
tion condition |[BJ, 



E 00 K, « E oq t e + eA max B = nh/2. 



(11) 
is the 



The range of directed areas — A max < A < A max 
product of the area L 2 of the billiard and the maximum 
number of times n ma x ~ vpT/L that a trajectory can 
encircle that area (clockwise or counterclockwise) in a 
time T. Hence A max = vpTL < vfTeL and we find 



TlT) 



E 00 (B) = E£ p n — -ev F LB. 



(12) 



We conclude that a magnetic field shifts the lowest 
adiabatic level downward by an amount ev F LB which is 
independent of te- Eq. (|12|l holds up to a field B : 



ad 



at 



which the lowest adiabatic level reaches the Fermi level, 



B 



ad 



nh 



nh 



2eA r 



2tecvfL 



(13) 



We have added the label "ad" , because the true critical 
field at which the gap closes may be smaller due to non- 
adiabatic levels below E^o . For B = 0, the ground state is 
never an adiabatic state |H| . In the next section we study 
the effective RMT, in order to determine the contribution 
from non-adiabatic levels (return times T > te). 
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E. Density of states 

The pair of quantization conditions JSJ) determines the 
individual energy levels with T < te and \A\ < A max = 
vfTeL. For semiclassical systems with L/X F 3> 1 the 
level spacing 6 of the isolated billiard is so small that 
individual levels are not resolved and it suffices to know 
the smoothed (or ensemble averaged) density of states 
Pad(E). In view of Eq. © it is given by 

Pad (E) =N [ * dT [ dAP(T,A) 



nh(m+ 1/2) + eAB 



(14) 



o 

3 L 

-a 
c 

o r 

o 

— 
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lead with 
time delay 




chaotic cavity 



in terms of the joint distribution function P(T, A) of re- 
turn time T and directed area A. In the limit te — > oo 
this formula reduces to the Bohr-Sommcrfeld quantiza- 
tion rule of Ref. pj for B = and to the generalization 
of Ref. 20] for B ^ 0. The adiabatic density of states 
(|14fl vanishes for E < E^ p . Its high energy asymp- 
totics (meaning E > E^ p , but still E <C A) can be 
estimated using P{T,A) = P(A\T)P(T) with the con- 
ditional distribution P(A\T) (which will be discussed 
in the next section) and the return time distribution 
P{T) = exp {-T/t d )/t d . One gets 



lim p ad {E) = - 

E— »oo 



,—tb/td 



1 + ^ 
TD 



(15) 



The limit l|15[) is less than the value 2/8, which also con- 
tains the contribution from the non-adiabatic levels with 
T > te- 



III. EFFECTIVE RANDOM-MATRIX THEORY 

The adiabatic quantization applies only to the part 
of phase space in which the return time T is less than 
the Ehrenfest time te- To quantize the remainder, with 
T > te, we apply the effective random- matrix theory 
(RMT) of Ref. 5]. The existing formulation HQ3 does 
not yet include a magnetic held, so we begin by extending 
it to non-zero B. 



A. Effective cavity 

The effective RMT is based on the decomposition of 
the scattering matrix in the time domain into two parts, 



S(t) 



S c z(t) if t<T E 

Sg(t) ]£t>TE. 



(16) 



The classical, short-time part S c i(t) couples to N c i scat- 
tering channels of return time < te , which can be quan- 
tized adiabatically as explained in the previous section. 



FIG. 6: Pictorial representation of the effective RMT of an 
Andreev billiard. The part of phase space with long trajec- 
tories (return time > te) is represented by a chaotic cavity 
with level spacing 8 c s, connected to the superconductor via 
a fictitious ballistic lead with N a g channels. The lead intro- 
duces a channel-independent delay time te/2 and a channel- 
dependent phase shift (f>„, which is different from the distri- 
bution of phase shifts in a real lead. 



The remaining 



N q = N - N cl = Ne 



-te/to = 



= N ( 



cir 



(17) 



quantum channels, with return time > te, are quantized 
by RMT with effective T^-dependent parameters. 

To describe the effective RMT ensemble from which S q 
is drawn, we refer to the diagram of Fig. El following Ref. 
[To). A wave packet of return time t > te evolves along 
a classical trajectory for the initial te/2 and the hnal 
te/2 duration of its motion. This classical evolution is 
represented by a fictitious ballistic lead with delay time 
te/2, attached at one end to the superconductor. The 
transmission matrix of this lead is an iV e fT x N c f[ diagonal 
matrix of phase shifts exp [i$(_B)] (for transmission from 
left to right) and exp [i$(— B)] (for transmission from 
right to left). The ballistic lead is attached at the other 
end to a chaotic cavity having an iV ff x N e R scattering 
matrix So with RMT distribution. The entire scattering 
matrix S q (t) of the effective cavity plus ballistic lead is, 
in the time domain, 



S q (t)=e l ^- B ^So(t~TE,B)e 



i$(B) 



and in the energy domain, 

S q {E) = e lETE ' h e l ^- B) S {E,B)e l ^ B) . 



(18) 



(19) 



The level spacing 8 c g of the effective cavity is increased 
according to 



8 eS /8 = N/N eS = e 



te/to 



(20) 



to ensure that the mean dwell time 2ttTi / N c f{S f[ remains 
equal to td, independent of the Ehrenfest time. 
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For weak magnetic fields (such that the cyclotron ra- 
dius mvp/eB 3> L), the phase shifts &(B) are linear in 
B: 

~ $(0) + B$'(0) = $(0) + diag 2 . . . <^ eff ] . 

. (21) 

The phases 0„ are the channel dependent, magnetic field 
induced phase shifts of classical trajectories spending a 
time te/2 in a chaotic cavity. 

The conditional distribution of directed areas A for a 
given return time T is a truncated Gaussian j^jl 0] , 



P(A\T) oc cx P (-A 2 Mg)0(A I] 
A 2 , oc v F TL 3 , 



14), 



(22) 



with the unit step function. This implies that the 
distribution P(<f>) of phase shifts 4> = eAB/fi for T = 
/2 is given by 



P(</>) oc exp 




101), (23) 
(24) 



The constant c of order unity is determined by the billiard 
geometry and B denotes the critical magnetic field of 
the Andreev billiard when te —> 0. Up to numerical 
coefficients of order unity, one has 0] 



Bn 



V F T D 



(25) 



B. Density of states 

The energy spectrum of an Andreev billiard, for ener- 
gies well below the gap A of the bulk superconductor, 
is related to the scattering matrix by the determinantal 
equation 



Det [l + S(E)S* (-E)] =0. 



(26) 



Since S c i and S q couple to different channels, we may cal- 
culate separately the contribution to the spectrum from 
the effective cavity, governed by S q . We substitute the 
expression (|19|l for S q , to obtain 



Det 



1 + e 2lETE/h S (E, B)n(B)SZ {-E, B)Q*{B) 
n(B) = e *>iB)-i*(-B) = diag[e 2 *,e 2l< 



. e 2 ^N e ff 1 



= 0, 
(27) 

(28) 



In Ref. 0] the density of states was calculated from this 
equation for the case B = 0, when $7 = 1. We generalize 



the calculation to B ^ 0. The technicalities are very 
similar to those of Ref. [2^ • 

The scattering matrix Sp(E 1 B) of the open effective 
cavity can be represented by [2J, |25| 

S (E, B) = 1 - 2mW T [E - H {B) + inWW T ] ^ W, 

(29) 

in terms of the Hamiltonian Hq(B) of the closed effective 
cavity and a coupling matrix W. The dimension of Hq is 
M x M and the dimension of W is M x N e g. The matrix 
W T W has eigenvalues MS c ff/n 2 . The limit M — > oo at 
fixed level spacing i5 e fi is taken at the end of the calcu- 
lation. Substitution of Eq. (|29|l into the determinantal 
equation (|27fl gives a conventional eigenvalue equation 



Det [E - H cff (B)} = 0, 

H (B) 
-H$(B) 

WW T smu WQ(B)W T 
WQ*(B)W T WW T sinu 



H cS (B) = 

w 



cos it 



(30) 
(31) 

(32) 



We have abbreviated u = Ete/%- 

The Hamiltonian Ho (B) of the fictitious cavity has the 
Pandey-Mehta distribution j2(j , 



P(H) oc exp 



AM 5: 



M 



x [(BeHitf + b-tQmHij) 2 ] J. (33) 

The parameter b S [0, 1] measures the strength of the 
time-reversal symmetry breaking. It is related to the 
magnetic field by 0] 
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(B/B ) 2 



(34) 



The ensemble averaged density of states p e e(E) is ob- 
tained from the Green function, 



Pes(E) = 

= ((z 



oil' I 



(36) 



where the average (• • • ) is taken with the distribution 
C3). Using the results of Refs. [H|23| we obtain a self- 
consistency equation for the trace of the ensemble aver- 
aged Green function, 



G 



f G u 


Gl 2 > 




{ G 21 


G22 j 





TrQ n TrC/12 
Tr Q21 Tr Q22 



(37) 



The four blocks refer to the block decomposition (|31|l of 
the effective Hamiltonian. The self-consistency equation 
reads 
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Gn — G22, G12G21 — 1 + G11 
/ E 

= iV eff 



e 2 *Gii + G12 sinw 



= JVcff 



_^ 1 Gi2 + e 2»^ c 21 ] + cos u + G?u sin it ' 

e-^Gn + G 2 isinM 



iV Bff 



2E T \Bq) 2 ) G21 + ^ i [e" 2 *Gi2 +e 2 *G 2 i] +cosu + Gnsinu' 



with the Thouless energy Et = ^/2te>. 

From Eq. (|35|) we find the density of states 



(38) 
(39) 

(40) 



Pes(E) 



SeS 



-Im 



Gn + 



TE 



G11 



isinu(G 2 ie 2i ^ 



G 12 e" 



td cosu cosu + Gn sinu + \Gxie" 2i< ^i + ^G 2 ie 2i ^ 



3=1 



(41) 



Because N e s ^S> 1, we may replace in Eqs. (|38H41|1 the 
sum £\ f(4>j) by / #-P(0), with P(0) given by Eq. J^J. 
In the next section we will compare the density of states 
obtained from (|38H41(I with a fully quantum mechanical 
calculation. In this section we discuss the low and high 
energy asymptotics of the density of states. 

In the limit E -> 00, E < A we find from Eqs. (|38H4l)|) 
that G12 = G21 oc 1/E — > while Gn — > — i. Substitut- 
ing this into Eq. (|41|l we obtain the high energy limit, 



lim p e g(E) = 




This limit is larger than 2/5 e g because of the contribution 
from states in the lead, cf. Fig. EI Together with Eq. i|15|) 
we find that the total density of states, 

p{E) = p eS (E) + p ad (E), (43) 

tends to 2/8 for high energies, as it should be. 

At low energies the density of states p c s(E) obtained 
from the effective RMT vanishes for E < Ef?„. In the 
limit te td the lowest level in the effective cavity is 
determined by the fictitious lead with return time rg. 
This gives the same gap as for adiabatic quantization, 

= = ^ (I " 2</W ) * S~E ' eVFLB > 

cf. Eq. IQ. The two critical magnetic fields Bf and 
Bf coincide in this limit, 

Dcff _ R ad _ nh _ o / TpL 

(45) 



cf. Eq. (|13|l . In the opposite regime of small T£i we find 
a critical field of 



Bf =B \l- , if T£ « y/LT D /v F , (46) 

which is smaller than i3^ d so £? c = -Bc ff - I n the interme- 
diate regime ^Ltei/ve ^ t~e ^ 773, the critical field i? c 
is given by 

B c = min (Bf ,Bf) . (47) 

We do not have an analytical formula for Bf in this 
intermediate regime, but we will show in the next section 
that Bf drops below Bf so that B c = Bf. 

IV. COMPARISON WITH QUANTUM 
MECHANICAL MODEL 

In this section we compare our quasiclassical theory 
with a quantum mechanical model of the Andreev bil- 
liard. The model we use is the Andreev kicked rotator 
introduced in Ref. . We include the magnetic field into 
the model using the three-kick representation of Ref. [27| , 
to break time-reversal-symmetry at both the quantum 
mechanical and the classical level. The basic equations 
of the model are summarized in Appendix [S] 

In Fig. [7| we show the ensemble averaged density of 
states of the Andreev kicked rotator and we compare it 
with the theoretical result 1|43|) . The Ehrenfest time is 
given by 0, Q 

T E = X- 1 [In (N 2 /M) + 0(1)] , (48) 

with M the dimensionality of the Floquet matrix. We 
neglect the correction term of order unity. The mean 
dwell time is Td = (M/N)tq and the level spacing is 
8 = (2ir I 'M)h I 'to , with tq the stroboscopic time. The 
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FIG. 7: (Color online) Ensemble averaged density of states p(E) of the Andreev kicked rotator. The dark (red) curves show 
the numerical results from the fully quantum mechanical model, while the light (green) curves are obtained from Eq. 1431 with 
input from the classical limit of the model. The energy is scaled by the Thouless energy Et = Ti/2td and the density is scaled 
by the level spacing 8 of the isolated billiard. The parameters of the kicked rotator are M = 2048, N = 204, q = 0.2, K = 200 
in panel a and M = 16384, N = 3246, q = 0.2, K = 14 in panels b, c, d. The three-peak structure indicated by the arrow in 
panels b, c, d is explained in Fig. |H] 



relation between B/Bq and the parameters of the kicked 
rotator is given by Eq. IjAlOjl . 

In Fig. U T£ < tjj and we recover the RMT result 
of Ref. [llj. The density of states is featureless with a 
shallow maximum just above the gap. In Figs. 0d, c, d 
te and td are comparable. Now the spectrum consists of 
both adiabatic levels (return time T < te) as well as ef- 
fective RMT levels (return time T > te). The adiabatic 
levels cluster in peaks, while the effective RMT forms the 
smooth background, with a pronounced bump above the 
gap. 

The peaks in the excitation spectrum of the Andreev 
kicked rotator appear because the return time T in Eq. 
jQ is a multiple of the stroboscopic time t 7J. The 
peaks are broadened by the magnetic field and they ac- 
quire side peaks, due to the structure of the area dis- 
tribution P(A\T) for T a small multiple of tq. This is 



illustrated in Fig. |S| for the central peak of Fig. [7| The 
distribution was calculated from the classical map (|A11|) 
associated with the quantum kicked rotator. The same 
map gave the coefficient c = 0.55 appearing in Eq. 1231) . 

In Fig. we have plotted the critical magnetic field B c 
at which the gap closes, as a function of the Ehrenfest 
time. For te <C to the Andreev kicked rotator gives 
a value for B c close to the prediction Bq of RMT, cf. 
Eq. (jAlOfl . With increasing te we find that B c decreases 
quite strongly. In the figure we also show the critical mag- 
netic fields B^ d for adiabatic levels and B° s for effective 
RMT. The former follows from Eqs. JTJJ and (CT^ . 



R ad_ 7r R /2td7o 



(49) 



and the latter from solving Eqs. l|38H40fl numerically. As 
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-1 1 

A ^ A max 

FIG. 8: Conditional distribution P(A\T) of directed areas A 
enclosed by classical trajectories with T — 2ro, for K = 14, 
q — 0.2 and to — 5to. The distribution was obtained from 
the classical map 1A111 at 7 = 0. Trajectories with T = 2ro 
give rise to a peak in the density of states centered around 
E/E T = (m + l/2)ivh/2T Q , cf. Eq. ifHJ. On the energy 
scale of Fig. |7| only the peak with m — can be seen, at 
E/Et = 2.5 7r ~ 7.9. In a magnetic field this peak broadens 
and it obtains the side peaks of P(A\2tq). 




0.2 0.4 0.6 0.8 1 1.2 

FIG. 9: Critical magnetic field B c of the Andreev kicked ro- 
tator as a function of the Ehrenfest time. The Ehrenfest 
time te = A -1 ln(N 2 /M) is changed by varying M and N 
while keeping q — 0.2 and td/to = M/N = 5 constant. For 
the closed circles the kicking strength K = 14, while for the 
squares from left to right K = 4000, 1000, 400, 200, 100, 50. 
The solid curve is the quasiclassical prediction l)47ft . The open 
circles are obtained from the closed circles by the transforma- 
tion \te —* Xte + 1.75, allowed by the terms of order unity 
in Eq. @SJ. 



already announced in the previous section, B^ d drops be- 
low Bf with increasing te , which means that the lowest 
level -Egap is an adiabatic level corresponding to a return 
time T < te- The critical magnetic field is the smallest 
value of B° s and B^ d , as indicated by the solid curve. 
The data of the Andreev kicked rotator follows the trend 
of the quasiclassical theory, although quite substantial 
discrepancies remain. Our quasiclassical theory seems to 
overestimate the lowest adiabatic level, which also causes 
deviations between theory and numerical data in the low 
energy behaviour of the density of states (cf. Fig. [7] pan- 
els c, d). Part of these discrepancies can be attributed to 
the correction term of order unity in Eq. I|48|l , as shown 
by the open circles in Fig. 

In the regime of fully broken time-reversal-symmetry 
the distribution of eigenvalues is determined by the La- 
guerre unitary ensemble of RMT [2^, |2^| . The ensem- 
ble averaged density of states vanishes quadratically near 
zero energy, according to 



P(E) 



1 - 



n(4:TrE/5)\ 
AnE/8 J 



(50) 



In Fig. ^]] we show the results for the Andreev kicked 
rotator in this regime and we find a good agreement with 
Eq. (fSTTfl for t e <C t d . We did not investigate the t e 
dependence in this regime. 




FIG. 10: (Color online) Ensemble averaged density of states 
of the Andreev kicked rotator for fully broken TRS. The his- 
togram shows the numerical results, while the curve is the 
theoretical prediction 150H of the Laguerre unitary ensemble. 
Both the energy and the density of states are scaled by the 
level spacing 5 of the isolated billiard. The parameters of the 
kicked rotator are M = 2048, N = 204, q = 0.2, while K was 
varied between 200 and 250 to obtain an ensemble average. 
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V. CONCLUSION 

We have calculated the excitation spectrum of an An- 
dreev billiard in a magnetic field, both using a quasi- 
classical and a fully quantum mechanical approach. The 
quasiclassical theory needs as input the classical distri- 
bution of times T between Andreev reflections and di- 
rected areas A enclosed in that time T. Times T smaller 
than the Ehrenfest time te are quantized via the adia- 
batic invariant and times T > te are quantized by an 
effective random-matrix theory with r^-dependent pa- 
rameters. This separation of phase space into two parts, 
introduced in Ref. [f| , has received much theoretical sup- 
port in the context of transport jl8l 1271 . 13CX l3lL 13^. 133 . l34| . 
The present work shows that it can be successfully used 
to describe the consequences of time-reversal symmetry 
breaking on the superconducting proximity effect. 

The adiabatically quantized and effective RMT spec- 
tra each have an excitation gap which closes at different 
magnetic fields. The critical magnetic field B c of the An- 
dreev billiard is the smallest of the two values B^ A and 
£?° ff . For relatively small Ehrenfest time te <C tjj the 
critical field _B° ff from effective RMT is smaller than the 
critical field B^ d of the adiabatic levels, so B c = B^r . 
This value B ^ is smaller than the value Bq of conven- 
tional RMT [ll|, because of the re-dependence of the 
parameters in effective RMT. For te ^S> td the two fields 
Bf and Bf coincide, but in an intermediate regime of 
comparable te and td the adiabatic value B^ d drops be- 
low the effective RMT value B° s . This is indeed what 
we have found in the specific model that we have inves- 
tigated, the Andreev kicked rotator 0- The lowest level 
has T < te for sufficiently large te and B. This is a novel 
feature of the Andreev billiard in a magnetic field: For 
unbroken time-reversal symmetry the lowest level always 
corresponds to longer trajectories T > te 0, and thus 
cannot be obtained by adiabatic quantization [5lll0|. 
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APPENDIX A: ANDREEV KICKED ROTATOR 
IN A MAGNETIC FIELD 

The Andreev kicked rotator in zero magnetic field was 
introduced in Ref. • Here we give the extension to non- 
zero magnetic field used in Sec. IfVI We start from the 
kicked rotator with broken time-reversal symmetry but 
without the superconductor. The kicked rotator provides 
a stroboscopic description of scattering inside a quantum 



dot. The propagation of a state from time t to time t + TQ 
is given by the M x M unitary Floquet operator F with 
matrix elements [27J 

F mn = (XUY*nYnX) mn . (Al) 

The three matrices X,Y, and II are defined by 

V - A „«( M 7/6ir) cos (2-n-m/M) ( \9\ 
1 mn — u mn^ > v-^^v 

Y - & p -i{M/X2n)V{2nm/M) (A^) 

n mn = A/- 1 / 2 e-" r / 4 exp [i(n/M)(m-n) 2 ] . (A4) 
The potential 

V(6) = K cos(7rg/2)cos(0) + — sin (7rg/2) sin (20) (A5) 

breaks the parity symmetry for q =/= 0. Time-reversal 
symmetry is broken by the parameter 7. For kicking 
strengths K > 7 the classical dynamics of the kicked 
rotator is chaotic. 

The Floquet operator (|A1|) describes electron excita- 
tions above the Fermi level. The hole excitations below 
the Fermi level are described by the Floquet operator F* . 
Electrons and holes are coupled by Andreev reflection at 
the superconductor. The NxM matrix P, with elements 

fnm - nm X I otherwise , > J 

projects onto the contact with the superconductor. The 
integer Lq indicates the location of the contact and N is 
its width, in units of Aj?/2. We will perform ensemble av- 
erages by varying Lq. The process of Andreev reflection 
is described by the 2M x 2M matrix 

-n (l-P T P -iP T P \ 

V = { _ ipTp j _ pTp J • (A7) 

The Floquet operator for the Andreev kicked rotator is 
constructed from the two matrices F and V 0, 

T = V 1 I' 1 [\ p^V 1/2 - (A8) 

The 2M x 2M unitary matrix T can be diagonalized effi- 
ciently using the Lanczos technique in combination with 
the fast-Fourier-transform algorithm 35] . The eigenval- 
ues e i£m define the quasi-energies e m £ [0, 2ir]. One gap 
is centered around e = and another gap around e = n. 
For N <C M the two gaps are decoupled and we can 
study the gap around e — by itself. 

The correspondence between the TRS-breaking pa- 
rameter 7 of the kicked rotator and the Pandey-Mehta 
parameter b for K ^> 1 is given by [27| 

Iim6 v / Mff = i • (A9) 

K^oo v 12?r K ' 
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Here Mjj is the size of the Pandey-Mehta Hamiltonian 
[26|. Comparison with Eq. (|34|) gives the relation be- 
tween 7 and the magnetic field B, 



M3/2 



7 



D 



(A10) 



In RMT the gap closes when B = Bo, so when 7 = 7o = 

37tM-V 2t o/^- 

For the quasiclassical theory we need the classical map 
associated with the Floquet operator (|A8|I . The classical 
phase space consists of the torus < 9 < 2ir, < p < 
6tt. The classical map is described by a set of equations 
that map initial coordinates (6, p) onto final coordinates 
(9',p') after one period To |27|. 

6 X = e±p/3-V'(e)/6-2ira e - L , 

Pl = pTlsm(0 1 )TV'(9)/2-6Tra pi , 

9 2 = 9 1 ±p 1 /3- 2ir<T 02 , 

P2 = Pi- 6na P2 , 

6' = 8 2 ±p 2 /3 + 1 sm(9 2 )/3-2ira g ,, 

p' = p2±lsm{6 2 )TV , (9')/2-6Tra' p . (All) 

The upper/lower signs correspond to electron/hole dy- 



namics and V'{&) = dV/d9. The integers ag and a p are 
the winding numbers of a trajectory on the torus. 

The directed area enclosed by a classical trajectory be- 
tween Andreev reflections can be calculated from the dif- 
ference in classical action between two trajectories re- 
lated by TRS, one with 7 = and one with infinitesimal 
7. To linear order in 7 the action difference AS* acquired 
after one period is given by [27| 



AS = 7 (cos 6*i — cos 82) 



(A12) 



The effective Planck constant of the kicked rotator is 
Tieff = 67r/M, so we may obtain the increment in directed 
area AA corresponding to AS from 

^BAA = = —7 (cos6>i — cos6> 2 ) . (A13) 
n n c ff 6n 

Since |cos#i — cos(9 2 | < 2, the maximum directed area 
A max acquired after T/tq periods is 



A — 9 



T h Pro" 



t eB V 2tz> 



(A14) 
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